Finite temperature study of bosons in a two dimensional optical lattice 
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We use quantum Monte Carlo (QMC) simulations to study the combined effects of harmonic 
confinement and temperature for bosons in a two dimensional optical lattice. The scale invariant, 
finite temperature, state diagram is presented for the Bose-Hubbard model in terms of experimental 
parameters - the particle number, confining potential and interaction strength. To distinguish the 
nature of the spatially separated superfluid, Mott Insulator and normal Bose liquid phases, we 
examine the local density, compressibility, superfiuid density and Green's function. In the annular 
superfiuid rings, as the width of the ring decreases, the long range superfiuid correlations start to 
deviate from an equivalent homogeneous 2D system. At zero temperature, the correlation decay is 
intermediate between ID and 2D, while at finite temperature, the decay is similar to that in ID at 
a much lower temperature. The calculations reveal shortcomings of the local density approximation 
(LDA) in describing superfiuid properties of trapped bosons. We also present the finite temperature 
phcise diagram for the homogeneous two dimensional Bose-Hubbard model. We compare our state 
diagram with the results of a recent experiment at NIST on a harmonically trapped 2D lattice [Phys. 
Rev. Lett. 105, 110401 (2010)], and identify a finite temperature effect in the experiment. 
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I. INTRODUCTION 

Much progress has been made in the last decade in 
the use of trapped, ultracold atoms for the optical lat- 
tice emulation of tight binding Hamiltonianai"— . In the 
bosonic case, Quantum Monte Carlo (QMC) simulations 
are possible on very large lattices (algorithms scale lin- 
early with the system size) and at low temperatures (no 
sign problem). As a result, close contact between ex- 
periments and theory has been possible, with success- 
ful quantitative comparisons of momentum distributions 
and phase transition critical points^"—. Recent QMC sim- 
ulations of the homogeneous system have refined early 
determinations^"— of the critical point for the ground 
state superfluid-Mott insulator boundary to very high ac- 
curacyii, and have also examined the finite temperature 
behavior at integer fillingii. QMC studies which include 
a confining potential and hence determine a "state di- 
agram" showing which phases coexist as a function of 
interaction strength and characteristic density have also 
been reported at low temperature^^ and compared to 
experimenlj^i^. 

In current experimentsp'"— li^"— ultracold atoms are 
first loaded in a harmonic trap, and then an external 
sinusoidal potential is slowly ramped up to create the 
optical lattice. By varying the depth of the lattice, the 
interaction strength {U /t) is changed, thereby driving the 
system through a quantum phase transition from a su- 
perfiuid (SF) phase to a Mott insulating (MI) phase. Al- 
though the initial system can be prepared at a relatively 
low temperature, the ensuing system after ramp- up of the 
lattice has a temperature which is usually higher due to 
adiabatic and other heating mechanism a ^ ' ^ . Recent 
experiments have reported temperatures on the order of 



T/t w 0.9/fcB^. At such temperatures, the effects of ex- 
cited states become important, motivating investigations 
into the finite temperature phase diagram, quantum ver- 
sus thermal transitions, the location of spatially sepa- 
rated phases in a trap at finite-T, etc. Aspects of finite- 
T effects on the MI have been discussed in Ref.— i^i^. 
Several QMC studies have also appeared at zero and fi- 
nite temperatures, dealing with quantum criticality and 
phase coexistence in a trap^"— . 

The goal in this paper is to provide a study of the com- 
bined effects of a confining potential and finite temper- 
ature on the state diagram of the Bose-Hubbard model 
in two dimensions, generalizing previous finitc-T QMC 
workii at fixed density, and ground state QMC studies 
in the presence of a confining potentia l^^'^^ . In addition 
to the trapped state diagram, we present the finite tem- 
perature phase diagram of the homogeneous 2D Bose- 
Hubbard model. Our results quantify the spatial inho- 
mogcncity/phase coexistence which is created by the trap 
and the interplay of thermal fiuctuations, which give rise 
to a third, normal (N) liquid phase in addition to the 
usual ground state SF and MI phases of the Bose Hub- 
bard model^^. We use an appropriately scaled measure of 
particle number, the 'characteristic density', which allows 
systems of different sizes to be compared^^, to present 
the scale invariant finite temperature state diagram. We 
compare our state diagrams with the results of a recent 
NIST experiment on a harmonically trapped 2D lattice^, 
and identify a finite temperature effect in the experimen- 
tal data. To better understand the nature of the trapped 
phases, we investigate the correlation function decay. In 
the annular superfluid rings, the correlation decay is dif- 
ferent from an equivalent homogeneous 2D superfluid - 
matching for a short distance and then falling off at a 
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faster rate for longer distances. At zero temperature, the 
correlation decay is intermediate between ID and 2D, 
while at finite temperature, the decay is similar to ID at 
a much lower temperature. In short, the width of the trap 
and the temperature determine the 2D to ID crossover. 
These findings provide evidence for the breakdown of the 
local density approximation (LDA) for the description of 
superfluid properties of trapped bosons. 

The article is organized as follows. In Sec. II, the Bose- 
Hubbard model^ and the observables used to monitor 
the state diagram are introduced. The QMC methodol- 
ogy is also briefly described. In Sec. Ill, we present the 
QMC homogeneous system phase diagram at finite-T. In 
Sec. IV, we analyze finite temperature effects in a har- 
monically trapped system, present the state diagram at 
zero and finite temperature, and compare them with a 
recent NIST experiment^. In Sec. V, we examine the de- 
pendence of spatial correlations in the trapped superfluid 
rings, and present evidence for the breakdown of LDA. 
Finally, we summarize our results in Sec. VI. 



II. MODEL AND COMPUTATIONAL 
METHODS 

Cold bosonic atoms in the lowest band of an optical 
lattice can be modeled by the Bosc-Hubbard model, 

+VTYrf n, ~ n^n, , (1) 

i i 

Here a|, Ui are the creation and annihilation operators 
of a boson at site i, rii = a|aj is the number of bosons, 
t is the strength of hopping^ between nearest neighbors 
(here on a square lattice), and U is the on-site repulsive 
interaction. Vr is the curvature of the external harmonic 
trap which introduces inhomogeneity to the system, rf = 
xf + yf is the distance from the center of the trap, Xi = 
dxi, where d is the period of the lattice, and i an integer. 
The hopping t, interaction U, and trap curvature Vt are 
tunable parameters in the experiment. 

In this paper we use the modified directed-loop 
algorithm^! for the QMC simulation based on the Feyn- 
man path integral^, since the simple application of 
directed-loop algorithm to the Bose- Hubbard model, that 
has the widest applicability, is not efficient especially for 
large U/t. 

To obtain the phase diagram, we examine three local 
observables - the density pi = ( n.i ) , compressibility Ki , 
superfluid density ps, and Green's function g{i,j)- To- 
gether, these variables are able to distinguish the spa- 
tially separated phases in a trap. The local compressibil- 
ity is a measure of the response of the local density to a 



change in the chemical potential, 

rrf^[K(TK(0))-(n,(T)}(n,(0))] (2) 

Here /i^ = /i — Vrr? is the site dependent local chemi- 
cal potential which decreases away from the trap center, 
and jS = 1/T is the inverse temperature. For a trapped 
system, a qualitative criterion for the appearance of a 
local MI is that the local density be an integer, as for 
the homogeneous case. At T = 0, a more precise crite- 
rion for identifying a local MI is to require that take 
on a value equal to the compressibility, k = dp/ dp, of a 
homogeneous system of commensurate fillingi^i^. It is 
typically the case that an extended (more than 4-5 sites) 
region of nearly constant, integer density has Ki close to 
its value in the homogeneous MI, thus providing a more 
rapid, visual means to identify a likely MI region. For 
finite-T, since the MI-N boundary is a crossover, there 
is no sharply defined phase boundary, and there is arbi- 
trariness in drawing such a boundary. In delineating the 
crossover boundary for the homogeneous and the trapped 
system at finite temperature, we have taken k < 0.02/C/ 
to identify the Mott region. A superfiuid region has a 
nonzero superfiuid density, while the normal region is 
identified as the region with zero superfluid density and 
incommensurate filling with k > 0.02/U. 

The Green's function, or one-particle density matrix, 

is useful in identifying SF regions. Specifically, a g{i,j) 
which falls off sufficiently slowly as |« — j| becomes large 
signals the development of off-diagonal long range order. 
Conversely, a rapid (exponential) fall-off oi g{i,j) is con- 
sistent with the appearance of a N or MI phase. 

To define a local SF density is more subtle, since ps in- 
volves the response of the system to a global phase twist. 
While there may be a way to generalize this to a local 
quantity, here we adopt the simpler LDA approach, as- 
signing the homogeneous system value of ps to different 
spatial locations of the confined systero^i. This gives an 
initial state diagram for confined systems. The correla- 
tion function g(i,j) provides more information about the 
superfiuid properties of the finite width annular super- 
fiuid rings. 

III. HOMOGENEOUS SYSTEMS AT FINITE 
TEMPERATURE 

The zero temperature ground state phase diagram 
for the homogeneous Bose-Hubbard model (Vr = 0) in 
two dimensions has been worked out using strong cou- 
plin g ex pansions^, and most recently using QMGii. In 
Ref. I III , a finite temperature phase diagram for a con- 
stant occupation p ~ 1 was presented showing the bound- 
aries of the SF and MI/N regions in the {T/t, U/t) plane. 
Mean-field theory has been uscd^"— to study different 
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FIG. 1: (color online). Occupation (p) and SF density (ps) as 
a function of the chemical potential {fJ./U) for two different 
temperatures T/t = 2 and T/t = 1/8 at t/U = 0.025 for the 
homogeneous 2D Bose-Hubbard model calculated using QMC 
simulations on a 32 x 32 lattice. A nonzero pa corresponds 
to a superfluid (SF) phase, zero ps with integer occupation 
p = 1 and vanishing compressibility to a Mott insulator (MI) 
phase, and zero ps with non-integer occupation to a normal 
(N) phase. The key point is that at the low temperature, 
wherever the density is incommensurate, the system is also 
SF. There are only MI and SF regions, and no N liquid. On 
the other hand, at higher temperature, normal regions develop 
which have incommensurate filling but also ps = 0. This is 
made clear in the figure at the incommensurate density region 
near p/U = where a rise in temperature from T/t = 1/8 
to T/t — 2 has destroyed the superfluid density turning the 
SF to N. Looking at the ps values, it is also evident that the 
superfluid properties at higher occupation (1 < p < 2) such as 
around p/U = 1 is stronger than that of the low occupation 
(0 < p < 1) around p/U = 0. 



properties and to obtain the finite-T phase diagram in 
any dimension. Various other approximate methods have 
also been used to study the finite temperature Bose- 
Hubbard model^i^. We present here the QMC homo- 
geneous finite-T phase diagram in the {pL/U,t/U) plane 
valid for the first two commensurate filling lobes. 

In Fig. [1] the density p and the SF density ps are 
shown as functions of p/U for two different tempera- 
tures T/t =1/8 and T/t = 2 for a 32 x 32 lattice at 
t/U = 0.025. (This lattice size is large enough so that fi- 
nite size effects are minor.) For T/t = 1/8, ps is non-zero 
whenever p is incommensurate. At this low temperature 
the system is either a MI or a SF. On the other hand, for 
T/t = 2, this is no longer the case. There are substantial 
regions where ps is zero even though p is incommensu- 
rate, which signifies a N liquid. For a homogeneous sys- 
tem, the boundary between the MI and normal regions 
in the finite temperature phase diagram is a crossover, 
as opposed to being a true phase transition at the SF- 
N boundary. The MI phase strictly speaking exists only 
at T = 0. However, Mott-likc features can persist at a 
finite temperature with a small value of compressibility. 
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FIG. 2: (color online). Finite temperature phase diagram 
for the homogeneous Bose-Hubbard model in two dimensions 
in the {p/U, t/U) plane for four different temperatures (a) 
T/t = 1/8, (b) T/t = 1, (c) T/t = 2 and (d) T/t = 4. 
(e) shows the phase diagram in a different representation, for 
constant T/U values of 0.02, 0.04 and 0.08. The lines that 
demarcate the SF and N is a phase boundary, and the lines 
that demarcate MI and N is a crossover. At finite tempera- 
ture, normal phase regions appear between MI and SF. These 
normal regions grow bigger for higher temperature. Although 
the finite-T homogeneous phase diagram is usually presented 
as in panel (e), with T in units of U , the phase diagram plots 
(a)-(d) for constant T/t values connect well to the confined 
system state diagram presented in Sec. IV, and provide a 
useful representation for optical lattice experiments. 



As mentioned in Sec. II, in drawing this MI-N crossover 
boundary at finite-T, we have taken k < 0.02/U to iden- 
tify the Mott region. 

The QMC homogeneous system phase diagram for the 
2D Bose-Hubbard model at finite temperature can be 
built up by sets of runs such as shown in Fig. [T] for differ- 
ent U/t, and is presented in Fig. [2j We show the phase 
boundaries for four different temperatures in panels (a)- 
(d); T/t = 1/8, which is a low enough temperature to be 
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the zero temperature phase diagram, T/t = 1, T/t = 2 
and T/t = 4. At T/t =1/8 the system has only two 
phases - MI with an energy gap inside the lobe, and a 
gapless SF outside the lobe. The transition between them 
is driven by quantum fluctuations as t/U is varied. The 
lowest value of U/t for which the SF-MI transition can 
occur is U/t K. 16.74 (i.e. t/U w 0.059)ii, the location 
being at the tip of the first lobe. At a finite tempera- 
ture, thermal fluctuations play a part, and MI lobes arc 
reduced. For example, at T/t = 1 in Fig. [IJb), we see 
that between the MI and SF phases, a region of normal 
phase has developed at the expense of reducing both the 
SF and MI regions. With further increase of tempera- 
ture, the normal phase region widens. For T/t = 2, the 
N-SF phase boundary around p = 1 is no longer lobe-like: 
the SF region ceases to bend back inward to low t/[/ at 
small density. With further increase of temperature to 
T/t = A, the SF between the first two MI lobes is also 
replaced by N liquid. 



In Fig. ^e) we present the finite-T phase diagram for 
temperatures in units of U, for T/U = 0.02,0.04 and 
0.08. As temperature increases, the MI lobes shrink, and 
the N liquid phase region grows between the MI and SF 
region. For T/U = 0.08 we find that the compressibility 
K > 0.02/(7 near t/U = 0, which is therefore no longer 
a MI in our criteria. The pictorial representations of 
how the MI lobes behave, although different in the two 
representations of T/t and T/U, are interchangeable. For 
a constant T/t phase diagram, as t/U decreases to zero, 
T/U also goes to zero, and the MI lobes in (a)-(d) for low 
t/U slowly approaches fi/U = 0, 1 and 2, giving rise to a 
lobe that looks pointy. Although finite-T homogeneous 
phase diagram is usually presented at constant T/U as 
in panel (e), phase diagram plots (a)-(d) for constant T/t 
values connect well to the confined system state diagram 
presented in Sec. IV, and provide a useful representation 
for optical lattice experiments where temperatures are 
often cited^— in units of T/t. 
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FIG. 3: (color online). The validity of the characteristic den- 
sity idea and the effects of temperature on spatially sepa- 
rated phases in a harmonically trapped system are illustrated. 
For T/t = 1/8, panel (a) shows the density profiles for two 
systems with the same characteristic density p = 34.5 but 
different trapping potentials and total number of particles, 
Vr/t = 0.025, Nh = 1380 and Vr/t = 0.05, Nt = 690. AU 
position dependent quantities are equivalenlii^ as a function 
of characteristic length scale ^ = ^t/Vr, as evident here in 
density matching, (b) shows the actual density profiles for the 
two trapping strengths, highlighting how these two different 
spatial profiles match in the scaled plot (a). Density profiles 
in (a) exhibit a constant integer plateau and an incommen- 
surate ring with n < 1, calculated using QMC for a trapped 
system, (c) shows the density profile at a higher temperature, 
T/t — 1, also showing a similar constant plateau and an in- 
commensurate ring. To determine whether the ring contains 
a SF or N liquid, we plot LDA derived SF density across the 
trap which shows that for (a) the entire ring is SF, whereas 
in (c) only a portion of the ring is SF and the rest is N: p is 
incommensurate but ps = 0. (d) shows the unsealed density 
profiles for the two trapping strengths in (c). 



IV. STATE DIAGRAM AT FINITE 
TEMPERATURE 

A. Characteristic density 



We conclude this section by noting that Fig. [5] can be 
reinterpreted as giving the local density and SF profiles 
in a confined system, assuming the validity of the LDA. 
That is, to each spatial site in a confined lattice are as- 
sociated the density and SF density of a homogeneous 
system with the same local chemical potential. The re- 
sulting sequence of LDA phases which occurs in moving 
away from the trap center can most easily be understood 
by starting with the location {n/U, t/U), where n is the 
chemical potential at the trap center, on the homoge- 
neous phase diagram of Fig. [5] and moving downwards, 
since decreasing fi/U corresponds to increasing Vt rf and 
thus moving outward from the trap center. 



For the Bose-Hubbard model in the homogeneous case 
(Vt = 0), the thermodynamic limit is reached by increas- 
ing the linear lattice size L while keeping the density 
p = N\)/L'^ constant, where A^f, is the particle number, 
and d is the dimensionality. The phase diagram is a func- 
tion only of U/t, T/t and p, not of Nh and L separately. 
In the same way, for the harmonically confined inhomo- 
geneous system with Vr 7^ 0, the generalization of this 
procedurei^ is to define a rescaled length, = Xi/(^ with 
^ = \Jt/VT, and a 'characteristic density', p = Nt,/^'^. 
For the 2D model, p = N^Vr/t. The state diagram and 
the properties of the inhomogeneous system then depend 
only on the combination p. Measurements at fixed U/t, 
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FIG. 4: (color online). State diagrams for a harmonically 
trapped system for (a) T = 0, (b) T/t = 1 and (c) T/t = 2. 
A schematic of all the phases according to their labels is shown 
in Fig. [S] The phase boundaries are drawn according to the 
multiple phases that coexist in a region. For example, in (a), 
region (II) corresponds to a state which has a MI plateau in 
the center surrounded by a SF shell. For higher T, normal 
(N) phase appears and therefore there are states with many 
varieties of coexisting phases. For example, phase (IV) in (b) 
has a MI in the center, surrounded outward by N, SF and N 
respectively. At a finite temperature, the tail always contains 
a normal liquid phase. For higher T, the SF-N phase bound- 
ary shifts to a lower value of U/t, and the first appearance of 
a MI plateau occurs at a higher value of U/t. 



T/t and p then approach a well-defined large L limit. Po- 
sition dependent quantities match when plotted in units 
of rescaled length {xi/£,). Any LDA-dcrived quantities 
like p or local ps trivially follow the rescaling. 



The validity of characteristic density idea discussed in 
Ref.lllis illustrated in FigM For T/t =1/8 (low enough 
to be effectively zero temperature), (a) shows the QMC 
density profiles for two systems with the same character- 
istic density p = 34.5 but different trapping potentials 
and total number of particles, Vr/t = 0.025, Ni, = 1380 
and Vr/t = 0.05, Nf, = 690. When plotted as a function 
of characteristic length scale ^ = ^/tJVr^ the density 
profiles are equivalent, as are all the position dependent 
local quantitiesi^. (b) shows the actual density profiles 
for the two trap strengths in (a), clearly showing the 
different extent of the profiles and highlighting the con- 
cept of equivalence in the scaled plots in (a). The same 
holds true at finite-T, as shown in Figs, ^c) and (d) for 
T/t = 1. The density profile in (a) shows a constant in- 
teger plateau and an incommensurate ring with ri < 1, 
calculated using QMC for a trapped system. In panel (c) 
the density profile at a higher temperature, T/t = 1, also 
shows a similar constant plateau and an incommensurate 
ring. To determine whether the ring contains a SF or N 
liquid phase, we plot LDA-dcrived SF density across the 
trap. This shows that for Fig. EJa) the entire incommen- 
surate ring is SF {ps ^ 0), whereas in Fig. ^c) only a 
small part of the ring is SF and the rest is N, where p 
is incommensurate but ps = 0. This is how we identify 
trapped phases containing distinct combinations of MI, 
SF and N regions. Due to the validity of characteris- 
tic density scaling, the state diagram we obtain below is 
scale-invariant - not depending on the specific values of 
trap strength (Vr), number of particles (Ni,) or lattice 
size, but on the combination NhVr/t = p. To obtain the 
state diagrams presented here, the largest lattice we used 
is 64 X 64. 



B. Zero temperature state diagram 

The state diagram for the harmonically confined sys- 
tem at zero temperature was presented in Ref. [l^ . In 
Fig. m^a), we show the T = state diagram in character- 
istic density - interaction strength (p, U/t) plane using 
QMC simulations in a trap, to higher values of p than 
presented in Ref. [l^. Three separate states arc possible 
here as indicated in Fig. EJa), and illustrated in Fig. \E\ 
The region (I) corresponds to SF phase all across the 
trap. In region (II), there is a MI plateau with (n = 1) 
in the center, surrounded by a SF ring of n < 1. (Ill) 
is a state (see Fig.[5jc)) where there is a local SF region 
at the center with n > 1, surrounded by a MI domain 
with n = 1 which is further surrounded by a SF ring 
with n < 1. This state diagram quantifies the parame- 
ter regimes for the appearance of these coexistent phases. 
Knowing the trap strength, total number of particles, and 
interaction strength we can predict the phase in a trap. 
A MI is obtained for a U/t that is always greater than the 
homogeneous system critical coupling of {U/t)c ~ 16.74. 
Only for a small window of p w 23, are the critical cou- 
plings comparable. A recent determination of the critical 
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FIG. 5: (color online). Schematic of all the phases for a trapped system at T = 0, T/t = 1, and T/t = 2, according to the labels 
shown in the state diagrams in Fig. [H At finite-T, part of the atomic cloud turns normal forming rings surrounding SF and 
MI regions, and three phases at T = increase to seven different phases at T/t = 1 and five at T/t — 2. The furthest part of 
the tail is always N at a finite-T as seen here in panels (d)-(o). Annular SF rings cannot form for a high enough temperature 
such as at T/t — 2 shown in panels (k)-(o). 



coupling by a NIST group^ can be understood in terms 
of the characteristic density trajectory the experiment 
follows when U /t ratio is increased^. Further experi- 
ments^ validate the agreement with the QMC trapped 
system state diagram. 

To go beyond the results of zero temperature state di- 
agram in Ref. HH, we compare T = QMC state diagram 
with a state diagram obtained by T = LDA. The green 
line in Fig. Hlja) is the LDA state diagram obtained from 
the homogeneous phase diagram by evaluating the den- 



sity profiles and characteristic densities for specific values 
of chemical potential at the phase boundaries^. The ob- 
vious disagreement is in the boundary between phase I 
(SF) and phase III (SF+MI+SF) which arises due to the 
finite gradients in the QMC simulations (and the exper- 
iments). Specifically, the disagreement between phases 
(I) and (III) boundary is related to the appearance of 
finite width MI shoulders. In the LDA picture the MI 
appears whenever we are at a chemical potential directly 
above the tip of the lobe in Fig. [21 On the other hand. 
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in a confined system, the shoulder develops more slowly 
as we move above {U /t)^ and only for a shoulder of suf- 
ficient width docs the confined MI compressibility equal 
that of a homogeneous system with the same U /t and in- 
teger occupation. The QMC boundary therefore slowly 
deviates from the straight line LDA boundary upwards. 
A recent 2D experiment done at NIST— has observed this 
deviation from the LDA. The transition from phase (II) 
to phase (III) requires the development of a superfluid 
bulge at the center of the trap. This occurs when it is fa- 
vorable to put a particle at the center rather than at the 
outer edge of the atomic cloud. By an energy matching 
condition, it is possible to show that p = ttU /<, which 
implies that the slope of boundary (II)- (III) should be 
TT. This is indeed the case for both the LDA and trapped 
QMC data points. For MI shoulders with n > 2 at higher 
/5, similar analytic phase boundaries exist with slopes Stt, 
9tt, etc. 



C. Finite temperature state diagram 

Figs, m^b) and (c) show the finite temperature state 
diagrams at T/f = 1 and T/t = 2 respectively. The 
phase boundaries demarcate the distinct states that are 
possible with different parts of the atomic cloud being a 
superfluid, normal fluid or Mott insulator. Fig. [5] illus- 
trates all the states of Figs. H] (b) and (c) according to 
their labels. In order to better understand all the phases, 
we can start at a low p at constant U/t = 40, and think 
what would happen if we put more particles in the trap, 
thereby increasing p along a vertical line at U/t ~ 40 in 
Fig. \Mjo). For low p < 1.35, N phase extends all across 
the trap (phase I). Putting more particles introduces a 
SF region at the center of the trap while the tail stays 
normal (phase II). Putting further particles in the trap 
introduces N liquid at the center while the SF-I-N regions 
of the previous phase remain, getting smaller (phase III). 
Adding more particles introduces Mott plateau at the 
center, with the N-f-SF+N region getting pushed out- 
side (phase IV). Following the phase schematics of Fig. [5] 
makes this process clear that adding more particles intro- 
duces a new phase region at the center while keeping the 
old phases, pushing them out, and reducing the extent of 
the regions. This process explains further appearance of 
phases V and VI with additional N and SF regions at the 
center. We can also think of varying U /t while keeping 
a constant p ~ 60. For U/t ~ 10, we start with a SF-|~N 
phase we encountered earlier. With increasing U/t we 
reach phase VII where the center is a SF surrounded by 
N+SF+N. 

With further increase of temperature to T/t = 2 the 
state diagram in Fig. \M,c) now has five phases as shown 
in Fig. [5l By focusing on a constant U/t ~ 40, and 
increasing characteristic density (by putting more parti- 
cles in a trap of fixed strength), we can visualize going 
through the phases, for this state diagram as well. First 
we encounter phase I with N liquid all across the trap. 



Further increase in number of particles introduces a MI 
at the center (phase II). Similarly, in phases III and IV, 
N and SF regions appear at the center pushing out the 
existing regions. If we start with N phase all across the 
trap (phase I) at a lower interaction U/t = 20 instead of 
U/t ~ 40, increasing p would introduce a SF bulge at the 
center (phase V). Note that T/t = 2 is high enough in 
temperature so that no confined SF rings can appear in 
the trap, as is the case for T = and T/t = 1. 

Some observations of finite temperature effects in a 
harmonically trapped lattice system are summarized be- 
low: (A) Thermal fluctuations always introduce a N 
phase region in the furthest tail of the atomic cloud where 
the density is small. As temperature increases, there al- 
ways appear N rings surrounding both MI and SF regions 
separately, thus giving rise to more than three conflned 
phases - seven at T/t = 1 and flve at T/t = 2. Thermal 
effects prevent SF rings to form for high enough temper- 
ature such as wc have seen in going from T/t = 1 to 
T/t = 2. (B) Finite temperature causes a shift of the 
SF-N transition coupling {U/t)c to a lower value than 
for a zero temperature trapped system, but to a higher 
value for the first appearance of a MI plateau. To see 
this, let us examine the zero temperature state diagram 
in Fig. Sl^a) where the phase boundary between states I 
and III separates the appearance of MI shoulders. As we 
raise the temperature to T /t ~ 1 such as in Fig. HKb), 
this boundary (SF-MI) moves to bigger U /t between the 
states VII and VI. While another boundary (SF-N) ap- 
pears between states II and VII at lower U/t that signi- 
fies the appearance of N phase in the center of the cloud. 
For temperature range T = to T/t < 1 as in NIST 
experiment^, the separation between these two bound- 
aries is small, and the MI transition would pass through 
a narrow range in U/t going through a SF-N transition 
first. In a {T/t ~ U/t) phase diagram at constant in- 
teger [n = 1), this lowering of U/t for SF-N transition 
can be understood in terms of a downward shift of tem- 
perature as has been observed and explained in a recent 
experiment^ in 3D. (C) In this paper, we identify and 
classify the trapped phases with QMC simulations, and 
do not study the issue of their experimental signatures. 
However, we would like to make some comments in re- 
lation to the state diagrams. There are two main ways 
of determining critical coupling in experiments - by an- 
alyzing momentum distribution (condensate fraction or 
peak width) in time-of-fiight (TOF) experiments^ and 
by direct (in situ) imaging of density profllea^^. For zero 
temperature trapped system, the critical coupling {U/t)c 
for SF-VII transition correspond to the appearance of a 
MI shoulder observed through a decrease in condensate 
fraction. At finite-T, the difference between coherence 
and incoherence as measured through condensate frac- 
tion differentiates between SF and N/MI. The distinction 
between N and MI may be difficult to detect with TOF 
imaging. For identifying local phases through the mea- 
surement of a global quantity, a larger portion (than the 
mere presence of N/VII) of the atoms has to turn N/MI 
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to detect the signature of its initial appearance. Since at 
a finite-T in a harmonic trap, there is a part of the cloud 
(in the tails) that is always normal and part of the cloud 
(annular rings) that is not a 2D superfluid as discussed in 
section V, signatures of the SF-N transition critical cou- 
pling (U /t)c in the momentum distribution of a harmoni- 
cally trapped system need to be carefully studied. In situ 
imaging may be a better method at finite-T to detect the 
appearance of normal phase and MI plateau. Single site 
resolution microscope a^^i^^ is a big step forward in this 
direction; detection of coexisting phases through density 
scanning has been discussed in Ref [24l 



In a recent experiment at NIST— , the state diagram for 
a 2D harmonically trapped lattice system was obtained 
by measuring the condensate fraction to identify the su- 
perfluid and Mott insulator regions as a function of atom 
density and lattice depth. By comparing to the T = 
QMC state diagram, they have shown the breakdown 
of LDA, indicated by the appearance of MI shoulders 
for higher values of {U/t)c than a homogeneous system; 
a phenomenon illustrated in our theory study here in 
Fig. Ha). 



D. Comparison with NIST experiment 
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FIG. 6: (color online). Comparison of theoretical state dia- 
gram boundaries with NIST experimental results of Ref. @. 
Circles and triangles represent respectively MI and SF data 
points from the NIST experiment; same set of data is plotted 
in (a) and (b). Comparison is made with theoretical state 
diagram boundaries in (a) T = and (b) T/t = 1. The types 
of coexisting phases for (a) and (b) are given in Fig. 3] There 
are MI points for low values of p in (a) that are below the 
state diagram boundary p « 18. For finite-T, the state dia- 
gram boundary goes down as in (b) containing many of these 
MI/N data points. We identify this feature as a finite temper- 
ature effect fully explainable with a finite-T state diagram. 



In Figs. |ni[a) and (b) we show respectively the T = 
and T/t = 1 state diagrams of Figs. |4|a) and (b) over- 
laid with the experimental data points of the NIST ex- 
periment; the circles represent MI and the triangles SF. 
We do not show the experimental phase boundary and 
its uncertainty reported in Ref. d. Although the experi- 
mental comparison was done with a T = state diagram, 
the temperature was reported to be T/t = 0.9(2), con- 
sidered a low enough temperature for such a comparison. 
However, as we have seen in Figs. Hl^a) and (b), because 
of the appearance of normal phase regions in the cloud, 
the state diagram at T/t = 1 begin to differ from T = 0. 
Here we choose to overlay the experimental data both in 
our T = and T/t = 1 state diagrams (same data in 
both panels in Fig. [5]) . It is visually quite clear that the 
experiment obtains SF points for low U/t and MI points 
for high U /t. However, the transition occurs at different 
values oi U/t for different p in the vertical axis, shown 
more clearly in the phase boundary plotted in Ref. [^, 
signifying breakdown of LDA. 



We find that for MI points in low p, the data agrees 
more with the state diagrams at T/t = 1 than at T = 0. 
For low p the measurement of different phases is diffi- 
cult due to increased thermal fluctuations, and the mea- 
surement uncertainty is large. In Fig. [Slla), there ap- 
pear many MI points below the boundary of p « 18. In 
Fig.EJb), we see that the lower boundary for MI/N goes 
down to p w 9 containing many of these points. We 
should emphasize that the reported experiment temper- 
ature of T/t = 0.9 is the initial temperature, while the 
final temperature during the measurement is unknown 
and could be higher, and in such a case, the boundary 
would go further down containing all the experimental 
MI/N data points. With the results of MI points below 
the lower boundary for T = 0, we would argue that the 
experiment was able to detect a finite temperature ef- 
fect, fully explainable with a finite-T state diagram. We 
would like to point out that with time of flight imaging 
of momentum distribution and measuring of condensate 
fraction, the distinction that is made in the experiment 
between SF and MI at T = is the same as the distinc- 
tion between SF and N/MI at finite temperature. 



9 




(b) 


* homogeneous 2D correlation decay 


• 2D trapped correlation at i=10 


t 









°0 5 10 15 20 25 



d 

FIG. 7: (color online), (a) Density profile and local super- 
fluid density for trapped atoms when the whole region is SF 
for parameters t/U = 0.025, T/t = 1/8, Vr/t = 0.012 and 
fJ'center/U = 0.1. In (b) clrcles show the correlation decay 
between a point i = 10 and other points along the circumfer- 
ence of radius i = 10, as a function of the arc distance d. The 
pentagrams show correlation decay for a 2D uniform system 
for a chemical potential same as that in the trapped system at 
i = 10. Since they match very well, we can conclude that the 
effect of trapping has not changed the long range correlations 
in this special case when the SF region is wide enough. 



E. Role of LDA-derived superfluid density in the 
state diagram 

The phase boundaries for finite temperature state di- 
agrams in Fig. m^b) and (c) were obtained determining 
the density profile and compressibility with QMC sim- 
ulations of a confined system. In addition to these two 
quantities, and in the absence of a definition of the local 
superfluid, LDA-derived ps was used to determine the 
trapped superfluid region as has been assumed in the lit- 
erature^iii^. In next section wc study bosonic Green's 
function which measures the off-diagonal long range or- 
der, and by comparing trapped correlations with homo- 
geneous (LDA) correlations, wc find that LDA does not 
correctly capture the trapped superfluid properties. For 
SF rings of narrow width, correlation decay do not match 
a 2D decay, and therefore assigning local ps values from 
a 2D homogeneous system has limited validity. So the 
state diagram boundaries we have obtained here give the 



upper limit boundaries of confined phases. This is only 
true for the regions where there are SF rings. High tem- 
perature state diagram such as at T/t = 2 in Fig. HJc), 
where no SF rings can form, the phase boundaries are not 
influenced by this. In Figs.|4l[a) and (b), for regions with 
SF rings, the superfluid designation with LDA-derived ps 
has to be kept in mind. At T = and T/t = 1, the SF 
rings are ID superfluids or a crossover between ID and 
2D. 

To be clear, there can be two LDA issues; one is due to 
finite size, and another one is in the description of super- 
fluid properties in trapped rings with LDA p^, an effect 
we identify in next section. Since QMC simulations in a 
trap were used for density and compressibility, we do not 
have the finite-size LDA issue related to the appearance 
of MI shoulders like in T = diagram in Fig. lUJa). It is 
the use of local LDA ps to designate the phases, for which 
wc would like to qualify our state diagrams by the fact 
that in the regions where there are annular SF rings, they 
can exhibit true long-range order or quasi-long-range or- 
der depending on the width of the ring. 



V. CORRELATION FUNCTION AT ZERO AND 
FINITE TEMPERATURE 

To further analyze the trapped system, we study spa- 
tial correlations in the trap by looking at the one-particle 
density matrix between bosons on site i and j, g{i,j) = 
{a\aj). As shown in previous section, the harmonic trap 
at zero and low enough temperature creates a superfluid 
ring surrounding the Mott insulating or normal region. 
The SF ring, having a finite width and a relatively longer 
length around a circle, has a quasi-lD geometry. 

For the parameters p./U = 0.1, t/U = 0.025, Vr/t = 
0.012, and zero temperature {T/t = 1/8), Fig.[7l^a) shows 
the density profile of trapped atoms where the whole re- 
gion is SF. The red circles in (a) are local superfluid den- 
sities obtained by LDA. In Fig.[7l^b), we plot the trapped 
correlation function decay (circles) between a point at 
J = 10 and all points along the ring at the same ra- 
dius, as a function of their distance d. Since the circular 
ring has the same local chemical potential at all points, 
we can take a homogeneous system with that specified 
chemical potential, and make a meaningful comparison 
with the trapped correlation function as we do in LDA 
for variables such as density, compressibility, etc. Such 
a homogeneous correlation decay shown in pentagram 
symbols match very well the trapped g{i,j) decay. This 
implies that the effect of trapping (inhomogeneity) has 
not changed the behavior of long-range correlations; this 
is true in this specific example where the SF region is 
wide enough that it retains its 2D superfluid properties, 
as we will soon show that trapping does have an effect 
on the SF rings of quasi-lD geometry. 

In Fig. m we investigate the zero temperature (T/t = 
1/8) correlations for a narrow SF ring. Panel (a) shows 
the density profile and local superfluid density indicating 
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FIG. 8: (color online). The effect of reduced dimensionality 
in a trapped superfluid ring. Panel (a) depicts a narrow SF 
ring of only 4 sites in width, for parameters t/U — 0.025, 
T/t = 1/8, Vr/t = 0.07 and Hcenter/U = 0.7. 3D plot of 
g{i,j) for i = 20 is shown in (b), where we see that g{i,j) 
is nonzero all across the ring. However, a comparison with 
the homogeneous correlations in (c) reveals that the trapped 
decay (circles) is faster in the ring than a true 2D decay (pen- 
tagrams), matching only for a very small distance. In (c) we 
also compare trapped ring correlations with a ID correlation 
decay (triangles) at the same parameters and temperature 
T/t = 1/8, and find that the ID algebraic SF decay is faster. 
This demonstrates that the bose gas in the incommensurate 
ring is in a crossover regime between ID and 2D. Assigning 
2D LDA ps values to the trapped rings therefore has limited 
validity, pointing out a breakdown of LDA in describing su- 
perfluid properties in a trap. 



that the ring is approximately 4 sites in width. Panel 
(b) shows a 3D plot of g{i,j) for a point in the SF ring, 
between i = 20 and for all other sites j in the lattice. It 



is evident that order persists all across the narrow ring 
as g{i,j) is nonzero along the circle, whereas it decays 
rapidly to zero radially from the SF region to the MI 
plateau. Similar to Fig.[71[b), in Fig.[5l^c) we compare the 
trapped correlation decay for i = 20 along the ring (cir- 
cles) to the homogeneous correlations at the same /i (pen- 
tagrams). It shows that the trapped SF decay matches 
the homogeneous decay for a short distance, after which 
it continues to deviate from a 2D SF decay. Thus the 
SF ring does not have 2D superfluid properties, rather 
it exhibits quasi-long-ranged correlations, influenced by 
the width of the ring. So the trapped annular SF decay 
not matching a 2D SF decay, we can find out whether it 
might match a ID correlation decay, since the ring has a 
finite width. In Fig. \Elc) we compare 2D trapped corre- 
lations along the ring with ID homogeneous correlation 
decay (triangles) at the same density and temperature, 
T/t = 1/8. In ID, the decay is algebraic (|i — p is 
the decay exponent, a positive real number) for super- 
fluid and exponential (e"'*"-''/^, ^ is correlation length) 
for a Mott or normal phased""—. We see in (c) that the 
2D ring correlations decay slower than a homogeneous 
ID algebraic SF decay but faster than a homogeneous 
2D decay. We can therefore conclude that the bose gas 
in the SF ring is in a crossover regime between a ID and 
2D superfluid. As the width narrows, it would approach 
a ID SF decay, and as the width increases, as we have 
seen in the earlier example, the trapped correlation decay 
would approach a 2D SF decay. 

While all these may seem quite intuitive because of the 
quasi-lD nature of the ring, the significance of these re- 
sults is in the fact that assigning 2D LDA-dcrivcd local 
superfluid density to the ring, as have been assumed in 



many recent paperi 



.24.28 



and in previous section, is not 



totally justifled. It indicates a breakdown of LDA for 
local condensate properties. We know that for trapped 
atomic systems, LDA describes very well the local den- 
sity dependent quantities, such as density, variance and 
compressibility^li^. A known failure of LDA in trapped 
systems is due to finite size effects in determining the 
appearance of a phase boundary^ii^. We show here that 
the effect of reduced dimensionality in traps lead to an- 
other failure - in the LDA description of local superfluid 
density in the trapped annular rings. 

Fig. HI illustrates the behavior of spatial correlations 
at finite temperature which is one of the key results of 
this article. For the parameters T/t = 1, n/U = 0.4, 
t/U = 0.05 and Vr/t = 0.4, Fig. E^a) shows the density 
profile and local superfiuid density exhibiting a super- 
fiuid ring of finite width. Fig. IHlJb) plots g{i,j) between 
the point i = 20 and other points on the lattice, showing 
clearly that although i = 20 is SF according to LDA pic- 
ture, the decay is such that there is no long-range order 
across the entire ring. Fig. Wic) quantifies the decay in 
greater depth - the circles showing the correlation along 
the ring slowly going to zero at a large distance. For a 
homogeneous 2D superfiuid with the same t/U , T/t, and 
at the same ^ as the spatial location i = 20, the correla- 
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FIG. 9: (color online). The behavior of spatial correlations for trapped superfiuids at finite temperature. For the parameters 
T /t = 1, n/U = 0.4, t/U = 0.05 and Vr/t = 0.4, (a) shows the density profile and local superfluid density indicating a superfluid 
ring of finite width, (b) plots g{i,j) between the point i = 20 and the rest of the lattice, showing clearly that although i — 20 
is SF according to LDA picture, the decay is such that there is no long-range order along the entire ring, (c) quantifies the 
decay in greater depth - the circles showing the correlation along the ring slowly going to zero at a large distance. g{i,j) for a 
homogeneous 2D superfluid with the same parameters is given in pentagrams. The triangles, ID correlation at the same density 
and temperature T/t = 1, decays much faster. However, the diamonds depicting ID correlations at a much lower temperature 
{T/t — 1/8) are seen to be similar to that of the 2D trapped annular SF correlations at T/t = 1. This suggests that for the 
purpose of long range properties the quasi- ID ring is at a much lower effective temperature than the equilibrium temperature 
of the 2D atomic cloud. For normal region at i = 15, the 2D homogeneous correlations match as shown in (d). 



tion is shown in pentagrams, where the decay is correctly 
that of a 2D SF. ID correlation at the same temperature 
T/t = 1 and parameters, shown in triangles, decays much 
faster. However, the diamonds depicting ID correlations 
at a much lower temperature T/t = 1/8 are seen to be 
qualitatively similar to that of 2D trapped annular ring 
at T/t = 1. This implies that as far as the long range 
properties are concerned, the annular ring is at a lower ef- 
fective temperature than the 2D cloud. This may be due 
to reduced fluctuations induced by the collective behavior 
of a finite width ring. In Fig. EJd) we show correlations 
in trapped normal region at i = 15 where it matches with 
the 2D homogeneous correlations. 



VI. CONCLUSION 

In this paper, we determined the combined effects of 
harmonic trapping and temperature for a square lattice 



Bose-Hubbard model to obtain the finite temperature 
state diagram. This extends previous worki^ which ex- 
amined phase coexistence at T = 0. As temperature 
increases, thermal fiuctuations melt away both the SF 
and MI phases, introducing the N phase. At finite-T, 
the N liquid phase is always present in a trap, in the 
lower density regions furthest from the trap center. Fur- 
thermore, each SF and MI region is surrounded by a N 
ring. This gives rise to many different confined phases. 
As the temperature increases, the critical coupling ([/ /t)c 
for the SF-N transition is lowered and the MI-N crossover 
coupling is increased. For the phases that contain SF an- 
nular rings in the state diagram, the quasi-long range 
nature of their correlations have to be kept in mind. In 
addition to the trapped state diagram, we presented the 
homogenous system phase diagram at finite temperature, 
for temperature in units of both T/t and T/U. 

We compare our state diagrams to a recent experiment 
at NIST— done on a harmonically trapped 2D lattice. 
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Although the experiment mainly focused on reporting 
the breakdown of LDA in obtaining the critical coupling, 
we identify that they have also observed signatures of 
finite-T in the state diagram. 

To further understand the trapped phases, we examine 
the dependence of spatial correlations g{i,j) in the annu- 
lar superfluid ring. We show that the correlation decay 
in SF rings docs not match the 2D homogeneous super- 
fluid at the same parameters. For short distances, on 
the order of the width of the ring, the correlations agree, 
while for longer distances the deviation gets bigger. At 
zero temperature, the correlation decay is intermediate 
between ID and 2D decay. At finite temperature, the 
trapped correlation decay rate is much faster than a ho- 
mogeneous 2D decay. Although it is still slower than a 
ID decay at the same temperature, ID correlations at 
a much lower temperature matches the trapped decay. 
This indicates that the ring is at a lower effective temper- 
ature, a fact that may have important consequences for 
long range properties in lower dimensions. These studies 
point out the fact that assigning LDA ps values to the 
trapped SF rings has limited validity, and thus provide 



evidence for the breakdown of the local density approxi- 
mation (LDA) in the description of superfluid properties 
of trapped bosons. 

The quantitative values for the phase boundaries pro- 
vided here provide numerical benchmarks for continuing 
efforts to emulate the Bose-Hubbard model on optical 
lattices, and demonstrate experimental-theoretical con- 
sistency for the numerical values of the location of the 
critical points. 
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